Multi-flavor bosonic Hubbard models in the first excited Bloch band of an optical 

lattice 
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We propose that by exciting ultra cold atoms from the zeroth to the first Bloch band in an 
optical lattice, novel multi-flavor bosonic Hubbard Hamiltonians can be realized in a new way. In 
these systems, each flavor hops in a separate direction and on-site exchange terms allow pairwise 
conversion between different flavors. Using band structure calculations, we determine the parameters 
entering these Hamiltonians and derive the mean field ground state phase diagram for two effective 
Hamiltonians (2D, two-flavors and 3D, three flavors). Further, we estimate the stability of atoms 
in the first band using second order perturbation theory and find lifetimes that can be considerably 
(10- 100 times) longer than the relevant time scale associated with inter-site hopping dynamics, 
suggesting that quasi-equilibrium can be achieved in these meta-stable states. 



I. INTRODUCTION 

The possibility to trap and manipulate the atoms in 
a Bose-Einstein Condensate using standing wave laser 
beams Q, |E 0, 0, IE El nas led to a renewal of the in- 
terest in basic solid state models. In such systems, the 
atoms experience a periodic potential from an optical lat- 
tice leading to formation of band structure in the energy 
spectrum. These bands have been investigated in exper- 
iments |7|. 

In the spectroscopy experiments in Ref. Q, the atoms 
experienced a periodic potential in only one direction, be- 
ing free to move on a much larger length scale in the other 
directions. This implied that interactions between atoms 
could be ignored. If the atoms are confined to reside 
on the sites of a lattice in three dimensions, interactions 
become important. As a result, it was shown theoreti- 
cally Q, and subsequently also experimentally [Bj, that 
a system of interacting cold atoms, residing in the low- 
est Bloch band of the periodic potential, maps onto a 
bosonic Hubbard model. This model is of great theo- 
retical interest since it exhibits a quantum phase transi- 
tion [E IE El El between ground states where the atoms 
are localized (Mott-Insulator) and where they are delo- 
calized (superfluid) as the strength of the hopping rel- 
ative to the inter-atomic interaction is varied. The dy- 
namics of particles under the influence of changes in the 
Hamiltonian (such as lattice tilts or rapid changes in the 
article interaction strength) has also proved interesting 

"innia. 

Another development is an interest in the idea of 
mixing bosonic atoms of different flavors in the lat- 
tice EE EE EE EE EEl • Several ways of achieving mul- 
tiple flavors have been suggested including using atoms of 
different species and exploiting different internal atomic 
states. 

So far, experiments on strongly interacting atoms in 
three dimensional optical lattices have been restricted 
to atoms in the lowest (zeroth) Bloch band. Recently 
Scarola and DasSarma considered the possibility of novel 



supersolid phases within the first excited Bloch band of 
an optical lattice. [2(j 

In this paper, the theory of atoms in the two lowest 
(zeroth and first) Bloch bands of a three dimensional 
optical lattice is considered. We show here, that due 
to the lack of available phase space for the decay prod- 
ucts, such excited states can (in some parameter ranges) 
have life times much longer than the characteristic time 
scales for inter-site hopping. Thus it should be possi- 
ble to establish quasi-equilibrium within the manifold of 
these metastable states. 

We find that it is possible in this way to realize novel ef- 
fective multi-species bosonic Hubbard Hamiltonians. De- 
pending on the choice of lattice depths the number of 
degenerate bands varies and we find effective models in- 
volving n flavors of bosons, where n can be I, 2, or 3. 
These flavors correspond to the three different possible 
nodal planes in the excited state wave function such as 
the one illustrated in Fig. . We will show that a char- 
acteristic of these Hamiltonians is that (to a good ap- 
proximation) each flavor can only hop in one direction 
(i.e., X (nodal plane) particles can only hop in the x 
direction, etc.). Neglecting interactions we would then 
have n interpenetrating one-dimensional free bose gases, 
one for each column (or row) in the lattice. Allowing 
intra-species interactions converts these one-dimensional 
gases into Luttinger liquids (or, if the interactions are 
strong enough, and the mean particle number per site is 
an integer, into Mott insulators). We show below that, 
besides intra-species interactions, the full interaction also 
includes on-site inter-species conversion terms that allow 
atoms to change flavor in pairs. Thus for example, two 
X particles constrained to move along a single x column 
can collide, turn into Y particles and move away along 
a y column. Such processes lead to novel quantum dy- 
namics for this coupled set of interpenetrating Luttinger 
liquids. 

As will be seen, the anisotropic nature of the hop- 
ping in conjunction with the pairwise conversion leads to 
Hamiltonians with an infinite but subextensive set of Z2- 
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FIG. 1: (Color online) On-site Wannier wave functions in the 
Harmonic oscillator approximation. The localized wave func- 
tions are to a good approximation described by harmonic os- 
cillator wave functions localized in each well. Above is drawn 
the wave functions 0(o,o,o) ( r ) (plotted in the plane z = 0) 
formed by the zeroth band Bloch functions and the wave func- 
tion 0(1,0,0) ( r ) formed by the zeroth band Bloch functions in 
the y- and z- directions and the Bloch functions from the 
first band in the x-direction. These are approximately har- 
monic oscillator states, 0(o,o,o)( r ) ~ exp[—a(x 2 + y 2 + z 2 )], 
and 0(1,0,0) ( r ) ~ x exp[— a(x 2 + y 2 + z 2 )] where the parameter 
a is determined by the curvature of the optical lattice poten- 
tial near its minima. Similarly 0(o,i,o) ( r ) ~ 3/0(0, o,o)( r ) an d 
0(o,o,i) (r) ~ 20(o,o,o) (r). 



gauge symmetries intermediate between local and global. 
Such infinite symmetries have been found in certain frus- 
trated spin models and in a 'bose metal' 
model j.25j and are known to cause dimensional reduc- 
tion in some cases. 0, m, 01 We will see below how 
this dimensional reduction appears in a simple way in 
this system. 

A related global Z2 symmetry and associated Ising or- 
der parameter appear in problems involving boson pair- 
ing due to attractive interactions mediated by Feshbach 
resonances. In that case the symmetry appears due to 
a conversion term that connects pairs of bosons with a 
distinct molecular field. This can lead to exotic states 
in which pairs of bosons are condensed but single bosons 
are not and in which half vortices are permitted prl |27| 

Further, due to strong interatomic repulsion, the 
ground state in 3D (three flavors) breaks a kind of chiral 
symmetry and displays an additional accidental ground 
state degeneracy at the mean field level. A similar sit- 
uation occurs for special parameter values in frustrated 
XY-models, where parallel zero energy domain walls can 
be inserted [2^]. The outline of this paper is as fol- 
lows: In section [H] the appropriate generalization of the 
bosonic Hubbard model is introduced along with numer- 



ical values of the parameters entering the Hamiltonians 
obtained from band-structure calculations for various lat- 
tice depths. Then, in section II I II the aforementioned 
effective Hamiltonians for atoms in the first band are de- 
rived for three particular choices of relative lattice depths 
in the xyz— directions. Using simple mean-field theory 
we sketch the ground state phase-diagrams in section llVl 
and in section[V]we discuss how the superfluid phases are 
reflected in the interference pattern in an experimental 
situation. Finally, in section IVT1 treating the interaction 
perturbatively to second order, we estimate the lifetime 
of a population inverted state (all atoms residing entirely 
in the first excited band). 



II. GENERAL LATTICE HAMILTONIAN 

The starting point is the Hamiltonian for weakly inter- 
acting bosons of mass m in an external potential |29| 



// = / d 3 x0t( x ) ( -^V 2 + Vo{*) + V r (x) ) $(x) 



1 4na s h 
2 m 



d 3 x0> t (x)^ t (x)V>(x)V>(x) 



(1) 



where a s is the s-wave scattering length. The external 
potential has two contributions Vo and Vr corresponding 
to the lattice potential and the magnetic trapping poten- 
tial. Denoting the wavelength of the lasers by A = 2a, a 
being the lattice spacing, the former can be written 



Vb(x)= F « si 



i—x,y.z 



Mir I '-^ x i ) . {%i}i=x,y,z = (X,y,z). 



The positional dependence of the magnetic trapping po- 
tential Vr — \m flfx" 2 , is much weaker than that 

i—x : y : z 

of the lattice, i.e. Hi <C x y ^^m 1 anc ^ wm ^ e ig nore d 
in the remainder of this paper. One should be aware 
though, that this term has been shown to influence for 
instance the phase diag ram of the single flavor bosonic 
Hubbard model j^Ef. 

For the cubic lattices considered here, the Wannier 
functions corresponding to the noninteracting part of the 
Hamiltonian in Eq. JQ) can be written 

n (x-R m )= ]~| <j>^J(xi-mid). 

i=x,y,z 

Here the bold face vectors n and m are integer triplets 
(n x ,n y ,n z ) and (m x ,my,m z ) which represent band in- 
dices and lattice sites respectively, i.e. 



m x ax + m y ay + m z az. 



These functions are to a good approximation described 
by localized harmonic oscillator wave functions sketched 
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in Fig. ^ The completeness of the Wannier functions 
allows the field operators to be expanded as 

V»(x) = E E ^n(m)</> n (x - R m ). 

m n 

The operators d„(m) and d n (m), which are the cre- 



ation and annihilation operators of bosons at site m 
and with band index n, obey Bose commutation rela- 
tions [d n (m),dl,(m')} = 6 n ,n'£m,m'- Ignoring all hop- 
ping other than nearest neighbor hopping and all inter- 
actions other than on-site interactions, the Hamiltonian 
in Eq. JQ| can be written 



6 w EE^H^H^H- E E*n E [4(m)d„(m')+4(m')dn(m) 

m n i=x,y,z n (m,m') j 

+ \ E E^( ni ' n2 ' n3 ' n4 ) Ki( m )4 2 ( m )^n 3 (m)d n4 (m) 



ni,n 2 ,n3,n4 m 



Here, the on-site interaction energies are defined as 



f/(ni,n 2 ,n3,n 4 ) 



At: a, % 



d d x 0* 1 (x)0* 2 (x)0 n3 (x)</>„ 4 (x) , 
I 



(2) 



(3) 



while the energies E n (m) and the hopping energies t„ 
are given by 

E n (m) ee J d 3 x^(x) (~|^V 2 + Vo(^j 0„(x) (4) 



where dimensionless Wannier wave functions 



*W EE / dx 



i)*( x .) { 



Voi(xi)) ^}{xi + a) 



(5) 

Note that the energies tn for hopping in the iEi-direction 
depend only on the lattice depth Voi in the corresponding 
direction and the i:th component ni of the band index n. 
The notation (m, m'). in Eq. J2J) indicates that the sum 
should be carried out over nearest neighbor sites m and 
m' in the a^-direction. One could for instance write, 



E =EE^.<w,^; + . 

) m m' 



It is straight forward to numerically solve the nonin- 
teracting Schrodingcr equation and find the energies in 
expressions ©-© above. In doing so, it is convenient 
to first switch to dimensionless units. Thus, we measure 
length in units of the inverse wave vector and potential 
depth in units of the recoil energy Er, i.e, £j = #i and 

V 0i EE Voi/ER, With E R EE £ (^) 2 . 

The hopping energies for the two lowest bands, ob- 
tained from band-structure calculations, are shown in 
Fig. |3 as functions of lattice depth. To get the on-site 
interaction (Eq. |3J) in a suitable form to use later on in 
the paper we define dimensionless overlap integrals 



o w (i;) = V27r / dei6»(«;0l IM«;OI' 



(6) 



have been introduced. The dependence on v in these 
functions is parametric, i.e. (f> n (vEn; £A/2tt) is the Wan- 
nier function corresponding to the one-dimensional non- 
interacting problem with potential depth vEr. The vari- 
able n denotes the band index. 

Approximating the Wannier functions with harmonic 
oscillator wave functions corresponding to the curvature 
at the bottom of each well, one finds approximate values 
for the overlap integrals 



,1/4 



(7) 



A comparison between these values for the overlap in- 
tegrals and those obtained from band-structure calcula- 
tions is shown in Fig. [3J 



III. EFFECTIVE HAMILTONIANS FOR ATOMS 
IN THE FIRST EXCITED BAND. 

In this section we will focus on the meta-stable sit- 
uation having all atoms in the first Bloch band(s) of 
the lattice. It is quite easy to achieve such a situation, 
i.e. consider an initial moment of time when the optical 
lattice has been loaded with atoms in the lowest Bloch 
band, n = (0,0,0). The anharmonicity of the lattice 
well potential allows one to treat the vibrational degree 
of freedom as a two level system. If one singles out, 
say, the ir-direction, then, by applying an appropriate 
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FIG. 2: Hopping energies and t^' in units of the recoil 
energy as functions of lattice depth Voi in the hopping direc- 
tion. The upper line is the hopping energy V[ for atoms in 
the first Bloch band hopping between nearest neighbor wells 
corresponds to atoms in the zeroth 



(i) 



while the lower line, 
band. 



,(t) 



FIG. 3: Overlap integrals O nn i(v) defined in Eq. ©. Solid 
lines, from top to bottom, Ooo,Ooi,On, obtained from nu- 
merical calculations. The dashed lines correspond to the val- 
ues in Eq. J7J obtained by using harmonic oscillator wave 
functions determined from the curvature of the potential at 
the well bottom. 



vibrational 7r-pulse, i.e. "shaking" the lattice in this di- 
rection with a frequency on resonance with the transition 
Tilj = £7(1,0,0) — ^(0,0,0) i the state can be inverted and the 
atoms excited to states with band index n = (1,0,0). A 
strong 7r-pulse can achieve this inversion in a time short 
compared to the inter-well hopping time so that the dis- 
persion of the upper band is not an issue in the inversion 
process. The simplest starting state would be the Mott 
insulator state with one boson per site in the lowest band. 
In a typical experimental setup, the parabolic confining 
potential will cause the population of each well to vary 
and the system will be in a state with regions of Mott 
insulators with different filling factors. It is however easy 
to confirm (by direct simulation) that even in this case, 
taking interactions into account, a pulse shape can be 
tailored that will invert the population simultaneously 
for regions with different filling factors provided a deep 
enough lattice is used. 

Another way of preparing the initial state is to use the 
method recently demonstrated by Browaeys et al. [3^ |. 
By loading a condensate into a moving ID lattice and 
applying a subsequent acceleration the condensate can 
be prepared in the lowest energy state (quasimomentum 
k = n/a) in the first Bloch band. The situation desired 
in this paper can then be obtained by ramping up the 
lattice in the two remaining (perpendicular) directions 
adiabatically. The natural question regarding the lifetime 



of the resulting meta-stable state will be considered in 
Sec. EH 

The subsequent dynamics of atoms in the first band(s) 
is then predominantly governed by some subset of the 
terms in Eq. J5J). This relevant subset will be referred to 
as the effective Hamiltonian. We will treat three different 
regimes of values for the lattice potentials Voj which lead 
to effective Hamiltonians with one, two, and three flavors 
respectively. The three scenarios are: 

1. Vox < V y,V Qz (ID, single flavor). 

2. V Qx = V 0y < V 0z (2D, two flavors). 

3- Vox = Vo y = Vqz (3D, three flavors). 

As indicated, the number of particle flavors as well as the 
dimensionalities in the effective Hamiltonians vary. 

The reason for the different numbers of flavors becomes 
clear if one considers the restrictions on the final states 
into which two atoms may scatter due to the interatomic 
interaction; The presence or absence of such states can 
be inferred from the presence or absence of degenerate, 
or nearly degenerate, levels in the energy spectrum of 
the noninteracting system. Take, for example the second 
scenario above with V$ x = Vq v <C Vq z and all atoms ini- 
tially in a state with index n = (1, 0, 0), then, due to the 
on-site inter-atomic interaction these atoms can scatter 
elastically into a state with index n = (0,1,0) through 
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a first order process the connecting different degenerate 
states. Further, it is easy to show that scattering result- 
ing in states with other indices, for instance n = (0, 0, 1), 
is only possible through higher order processes if energy 
(and also parity) is to be conserved and can safely be 
ignored if the gas is dilute. Hence, the atoms can, at 
a formal level, be divided into two flavors: an X flavor 
corresponding to atoms in n = (1, 0, 0) and a Y flavor in 
n = (0, 1, 0). By the same argument one can see how the 
one- and three-flavor situations arise. 

Apart from having different number of flavors the di- 
mensionalities of the effective Hamiltonians differ. To 
understand this consider again the second case above, 
Vox = Voy <C Vq z , with particles in the excited bands 
n = (1, 0, 0) and n = (0, 1, 0) corresponding to X- and Y- 
fiavors. For the X flavor, hopping in the x-direction has 
a matrix element {Vq x ) while hopping in the y- and 

z- direction have matrix elements t^\Vo x ) and tQ Z \Vo z ) 
respectively. Looking at Fig. [5] it is then clear that, to 
a good approximation, the X particles can only hop in 
the x-direction while hopping in the y- and z-directions 
is strongly (exponentially) suppressed. Similarly the Y 
particles can only hop in the y-direction and all hopping 
occurs only in the x — y plane, hence the 2D character. A 
similar argument holds for the three flavor case where in 
addition to the X and Y particles, there are Z- particles 
hopping in the z-direction. 

The effective Hamiltonians also contain terms arising 
from the on-site interaction. Apart from the terms that 
repel atoms from each other, the symmetry of the on-site 
interaction allows, say, two X particles moving in the x- 
direction to collide and convert into two Y particles which 
thereafter move off in the y-direction. The time reversed 
process can of course also occur. Thus, the number of 
particles of each flavor is not conserved and there is a 
pairwise exchange of particles of different flavors. 

The anisotropy of hopping and the flavor conversion 
process is schematically depicted for the 2D (two flavors) 
case in Fig. Particles of X flavor are shown in gray 
while the Y flavor is drawn in black. 

Below, we give the effective Hamiltonians for all three 
different cases listed above. 



A. ID Hamiltonian, single flavor, (Vox <C Vo y = Vq z ) 

The first case to be considered is when Vo x <C Vo y = 
Voz and only states with band index n = (1,0,0) (and 
possibly some residual atoms in n = (0,0,0)) are occu- 
pied. In anticipation of the other effective Hamiltonians 
it is convenient to introduce for the X flavor the creation 
and destruction operators X* and X, i.e.: 

X m = d(i, ,o)(m), = dj 10>0) (m) 

"m = 4),o,o) ( m Mo,o,o) (m). (8) 



(a) 



mi- 



(b) 



FIG. 4: (a) The Hamiltonian in Eq. li lt describes a 2D sys- 
tem where atoms, which can formally be thought of as having 
two different flavors (same type of atoms but in different local- 
ized on-site orbitals), hop around subject to on-site repulsive 
interactions. One flavor, the X-flavor, can only hop in the 
x-direction whereas the other, Y-flavor, can only hop in the 
y-direction. (b) Conversion process. The conversion term in 
Eq. 11 II takes two Jf-atoms on the same lattice site and turns 
them into two Y-atoms or vice versa. 



The effective Hamiltonian is then essentially that of a 
quasi one-dimensional bosonic Hubbard model 



Hid = 5^n&> [E x (m) + U 

- t 



U, 



E [xLx m , + h.c. 



(x) 

The energies t = t\ , Uq x , and U xx arise from the inter- 
well tunneling and the inter-atomic interaction respec- 
tively. The presence of atoms residing in the lowest band 
leads to an additional effective on-site energy and can be 
absorbed in the on-site energies E x (m), (s = x, y, z). Al- 
though this single-flavor model is equivalent to a single- 
flavor model in the zeroth Bloch band, the additional 
random on-site potential resulting from residual atoms 
could be exploited in the study of the disordered Bose- 
Hubbard system. 

The parameters entering the Hamiltonian © are con- 
veniently expressed as 

U xx = 2V2ttE r J Ooo(voy)0oo(vo2)Oii(voz), ( 9 ) 



U 0x = 2U xx O i (v 0x )/On(v 0x ) 



(10) 



B. 2D Hamiltonian, two flavors, (Vox = Vo y <S Vbz) 

To simplify the notation for the case Vq x = Voy *C Vq z , 
we introduce new letters for the creation/annihilation 
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operators Y m = d( ,i.o) ( m ); = Y^Y m . Then, the particles in any given column of the lattice running in the 
Hamiltonian governing atoms in the excited bands be- x-direction (and similarly for Y and Z particles). These 



comes 



U s . 



h™ = E E^m^+ E E« 



E [^m^m' + h.C.] - t J2 [ X L X ™' + h.C. 
(m,m') y 

+^/E ft 1 



(m,m') . 



(ii) 



Again, the energy U xy arises from the interatomic in- 
teraction and depends on the lattice depth. Note that 
this two flavor bosonic Hubbard Hamiltonian differs in 
an important aspect from previously studied two flavor 
systems: the presence of the last term that mixes the 
two flavors. Hence, the inter-atomic interaction leads to 
a "Josephson term" that allows for the conversion of two 
X-atoms into two Y-atoms and vice versa. The coeffi- 
cients U yy — U xx are given by the same expression as in 
the ID case while 

U xy = 2s/^E R (^) O o(voz)0oi(v 0x ) 2 . (12) 

Figure 0| illustrates the dynamics in the 2D (two flavors) 
situation. 



C. 3D Hamiltonian, three flavors, (Vox = Vo y = Vo z ) 

The generalization of the above Hamiltonian to the 
case when Vq x = Vo y = Vq z is straight forward. Intro- 
ducing a third flavor Z m = d( 0i0) i)(m), hm = Z^Zm, 
one may write an effective Hamiltonian as 



s—x,y,z m ^ 

J2 E Us »' ("m^ + ^[sLsLs'mS'm + h.C.] J 
s^s' m 



* E E [- 



h.c] 



(13) 



s=X,Y,Z <r 



Here V r ss > = 5 SS >U XX + (1 - 5 ss >)U xy with U xx and U xy 
given by © and 11211 with vq x = vq v = vq z . 



D. Z-z Gauge symmetry 

Because of overall number conservation the Hamilto- 
nian has the usual global C7 (1) symmetry. However, be- 
cause the flavor conversion occurs pairwise and locally 
(i.e., on site), the Hamiltonians described above also ex- 
hibit an infinite number of Z 2 gauge symmetries corre- 
sponding to conservation modulo 2 of the number of X 



symmetries correspond to invariance under each of the 
transformations 



U 



x 



U 



(m x ,m z ) 
Y 



exp 



exp 



ITT 



/ i (m x ,m B ,m; 



^E^i 



V 



U { "^ = exp 



^"E ^(mx,m ! ,,TOz)'^( m *' m !<> m z) 



where the integer pair (rrii, rrij) in the superscript of each 
U determine the location of a column. The first transfor- 
mation for example takes X^^^ -> -X^^^^j 
for all m x in the column specified by m y and m z . Since 

X' and X operators always appear pairwise, the Hamil- 
tonian is invariant under this class of Z 2 transformations. 
These Z 2 symmetries are in a sense intermediate between 
local and global. While the number of such symmetries 
is infinite (in the thermodynamic limit) it is of course sub 
extensive and thus not large enough to fully constrain the 
system (or to make it integrable for example). As men- 
tioned in the introduction, such symmetries have been 
found in certain frustrated spin models 0, |2^, US 0] 
and in a 'bose metal' model |25( and are known to cause 
dimensional reduction in some cases. [2ll I22I |24| Be- 
cause introducing a defect across which the sign of the 
Z 2 order parameter changes along any given single col- 
umn costs only finite energy, the system will, like the ID 
Ising model, disorder at any finite temperature thereby 
restoring the Z 2 symmetry. We will see below how this 
reduced dimensionality physics appears in a simple way 
in this system. 



IV. MEAN FIELD THEORY PHASE 
DIAGRAMS FOR THE EFFECTIVE 
HAMILTONIANS 

Having derived effective Hamiltonians in one, two, and 
three dimensions, we turn now to the investigation of 
their ground states. The ID, single flavor Hamiltonian 
has been extensively studied (see for instance Ref. |3^| 
and references therein) and needs no further discussion 
here. The other two Hamiltonians in Eqs. (|ll|l and i|13|) 
deserve some attention though. 



A. Phase diagram 2D, two flavors 

The 2D Hamiltonian i|ll|) is a two flavor bosonic Hub- 
bard Hamiltonian, a system that has recently received 
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much attention and shown to have a rich phase di- 
agram [3 E E E3 E H3- In this section we 
will investigate the ground state of the Hamiltonian in 
Eq. Ijlll) using simple mean field theory. The Hamilto- 
nian (|llf> differs from those previously studied in two 
aspects: the presence of pairwise inter flavor mixing and 
the anisotropic tunneling. 

We follow here the method suggested in Ref. [l^ (see 
also Refs. [!,E3)- We consider the possibility that the 
global U(l) and columnar ^-symmetries discussed in 
section IHI Dl are spontaneously broken by introducing 
complex scalar columnar order parameter fields ip x (rn y ) 
and ip y (m x ), i.e. one for each x-column and one for each 
y-column. These fields should then satisfy the self con- 
sistency conditions 

i>x(m y ) = (X( mxtmy )\ (14) 



for all m x in the x-column specified by m y , and, 



ip v (m x ) = {Y {m!C ,m v )) (15) 



for each m y in the y-column specified by m x . For sim- 
plicity we omit in this discussion of the 2D (two flavors) 
case the z component m z of the position vector m. The 
possibility that fluctuations restore the symmetry will be 
discussed further below. 

Mean field theory results from decoupling sites in the 
same column by neglecting fluctuations in the kinetic en- 
ergy. For instance, for the x-column specified by a par- 
ticular value of m v one has 



$tm*,m v ) ~ ^x( m y) + C( m y))(^(™ x +l,m y ) " lpx(m y ) + tp x (m y )) 
i>x(m y )xl m ^ my) +ip*(my)X {mx+1<mv) - \ip x (m y )\ 2 , 



Thus the sites along each column decouple. Doing the same for the Y:s and writing the Hamiltonian in dimensionless 
form where all energies are scaled by U xx , i.e., /i2«j = H2d/U xx , t = t/U xx and U xy = U xy /U xx , we obtain h,2D ~ 
S ^2d F ( m ' x l } x{fny)i ^("ij)). Here, the on-site mean field Hamiltonians are given by 

m 

h^/(m;ip x (m y ),^y(m x )) = -2t\^ x (m y )X^ + ijj* x {m y )X m - 2t i) y (m x )Y^ + ip*(m x )Y m +U xy n^n^ 



E - 1] - ^ + 2i|V^(m)| 2 ) + EnL [X lxlY m Y m + >,;,>,;, .V.„ .V 



(16) 



where fl = [m/U xx serves as a common chemical potential. 
The on site Hamiltonians satisfy the eigenvalue relations 

for two arbitrary complex fields. An eigenstate of the full 
mean field hamiltonian can be written as a product state 
of such eigenstates 

i*) = n i e « m (^( m ?/)'^( TO z))) 

Til 

where the fields satisfy the self consistency conditions in 
Eqs. (|14f> - l|15fl . The mean field ground state is obtained 
by globally minimizing the energy 

E = E e » m (Vk (m y ) , ip y (m x )) 

m 

with respect to the fields and the set of eigenstates {n m }. 
This is most easily done by numerical diagonalization in 



a truncated Hilbert space where each site can hold at 
most a total of -/V max atoms. Since 

min e nm (t/j x (m y ),'ifjy(m x )) 

[n m ,^x(m y ),V H (m x )] J 

> V" min e n (ip x ,4> y ) 

it is enough to minimize the ground state energy of a sin- 
gle site with respect to the fields and then find the largest 
manifold of states compatible with having columnar or- 
der parameters fields. Carrying out this scheme reveals 
two different scenarios for the minimum of each on-site 
energy e n (tp x ^ y ); either ^ x = tp y = and n x + n y is 
integer (incompressible), or, \ip x \ = \ip y \ ^ (compress- 
ible). The former case corresponds to a Mott insulating 
state while the latter suggests a superfluid phase. 

Due to the positivity of U xy the last term in the mean 
field Hamiltonian is minimal whenever ip x and ip y , on the 
same site, have a phase difference of ±7r/2. 
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For the (mean field) ground state manifold we must 
have in this phase \ip x (m y )\ — \ip y (m x )\ for all x— and 
y— column order parameters. Requiring the phases of all 
ip w in each x— and all ip y in the y— columns to be the 
same while fixing the relative phase between ip x and ip y to 
±7r/2 results in configurations as the one shown in Fig.0 
Here the phases of ip x and ip y are shown represented as 
arrows (planar spins) . The direction of the arrows defin- 
ing the angle. Clearly, this phase shows a breaking of the 
global U(l) symmetry. The meaning of the quasi- local 
nature of the Z 2 symmetries discussed above becomes 
clear. Although the phases of ip x in each x— column are 
the same there is no energy cost associated with flipping 
all the spins x-spins in a single column or all the y-spins 
in a y-column. The ordering between different columns 
is thus nematic. 

One should note here, that since the only energy cost 
associated with flipping a single spin, say an cc-spin in 
an column, is given by the states of the neighboring x 
spins in the column, the situation is essentially that of a 
ID Ising model along each column. Hence, at any finite 
temperature, domains of flipped spins will proliferate and 
the Z 2 symmetries will be restored. This essentially one 
dimensional behavior is an example of the dimensional 
reduction mentioned above. 

The model under consideration is highly anisotropic. 
Mainly since X particles can only hop in the x-direction, 
it seems to be impossible to develop phase coherence 
among X particles in different x-columns (and similarly 
for the other flavors). Suppose however that, as discussed 
above, the flavor exchange interaction term causes the 
relative phase of two flavors, say X and Y, to lock to- 



i 



4fy" 4jv 

4^- 4^» »J> 

4^» 4^» »J> 

<f- -f* 



FIG. 5: Columnar phase ordering in 2d superfluid phase. 
The directions of the arrows correspond to the phase angles 
4>x(m y ) and cf> y (m x ) of the order parameter fields ^( m j) = 
\iP x (m y )\e i4 '* < - m y ) and ip v {m x ) = \ip y (m x )\e i ' l '» ( - mm) . Solid ar- 
rows correspond to <f> x and dashed to <j> y . 



gether so that Y^X condenses 

V> = (y^x\ ¥= o. 

In this case the mean field decomposition of the exchange 
interaction yields terms of the form 

V ~ i)Y* It + tp* 

which permit individual particles to change flavor and 
hence phase coherence can freely propagate in all direc- 
tions throughout the lattice via a kind of 'Andreev' pro- 
cess (i.e. self-energy off-diagonal in flavor index) in which 
an X particle can turn into a Y particle when it needs 
to travel in the y direction. 

To understand this isotropic superfluid phase, it is 
convenient to consider a phase only representation with 
compact phase variables on each site X m — > e~ lv ™ and 
Y m — > e~ t(p ™. The flavor exchange ('Josephson') term 
then becomes (for the 2D (two flavors) case) 



V 



CW5^cos(2[v4 



Defining ip^ = ip^ ± ip^ we have 



V = U, 



x v 



cos(2vJ m ). 



Assuming that the relative phase of the condensates is 
locked together by this Josephson term is equivalent to 
assuming that (for U xy > 0) the fluctuations of ip~ away 
from the ground state value 7r/2 (or its equivalent — 7r/2 
under the Z 2 gauge symmetry) are massive and can be ig- 
nored. Thus we obtain tpm V ^ — Vm/^^ 71 "/^ The contin- 
uum limit of the anisotropic kinetic energy T = (d x ip x ) 2 + 
(d y ip y ) 2 then becomes T ~ ((d x ip + ) 2 + (d y ip + ) 2 ) and we 
immediately sec that the anisotropy has effectively dis- 
appeared at long wavelengths and we have a superfluid. 

A vortex configuration in ip + can be viewed as a bound 
state of two half vortices in the ip x and ip y fields. The 
columnar Z 2 symmetry allows the i\j x field to have a phase 
jump of 7T across a cut parallel to the x axis and similarly 
for Thus half vortices are permitted. If the two order 
parameter phases are locked together (<p~ fluctuations 
are massive) then the two half vortices are confined to 
each other as shown in Fig. ©. Such a vortex has an 
energy which scales (as usual) only logarithmically with 
system size, despite the semi-infinite branch cut (tt phase 
jump) of ip x running horizontally out to the right from 
the vortex center and of the similar branch cut in ipy 
running vertically out above the vortex center. To see 
that such a vortex is topologically well defined despite the 
Z 2 symmetry one can consider a loop around the vortex 
core as shown in Fig. ©. In going around the loop we 
add up the phase twist A(p and map onto the complex 
plane. To calculate A<fi along the loop the changes in (p v 
has to be added when going vertically and the changes 
in tp x when going horizontally. The net results is that 
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FIG. 6: (color online) Phase configuration for ip x (black solid 
arrows) and %p y (red dashed arrows) containing a half vortex. 
Notice the branch cuts indicated by the long solid black and 
dashed red lines. The Z2 symmetry means that these branch 
cuts have zero 'string tension' and contribute only a finite 
core energy to the vortex. The half winding number can be 
seen by going around a loop A-B-C-D-A and calculating the 
total phase twist A<f>. This twist is calculated by summing 
the changes in cp y when going vertically and ip x when going 
horizontally. 
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FIG. 7: Mean field ground state phase diagram for the 2D 
(two flavors) Hamiltonian in Eq. 1161 in the plane of fx, the 
scaled chemical potential [i/U xx , and t — t/U xx , the scaled 
hopping energy. Here calculated for the experimentally rele- 
vant ratio U X y = U xy /U xx = 1/3 using a truncated Hilbert 
space with at most 10 particles per site. Lobes of Mott- 
insulating states, of successively increasing integer filling fac- 
tor with increasing chemical potential jl, are surrounded by 
a superfluid phase. The superfluid phase is characterized by 
columnar order parameter fields ip x (m y ) and ip y (m x ), one for 
each x- and y- column respectively. All ip x and ip y have equal, 
nonzero, magnitudes while their relative phases are either 
or 7r. (see Fig. 



going once around the vortex core the phase winds by n. 
If one applies a 7r-flip in all the (p x (<p v ) phases in any 
row (column) the mapping onto the complex plane will 
remain invariant. 

The Mott insulating states, having integer number of 
atoms in each well, are best characterized by the t = 
eigenstates. These are product states 

|*o(i = 0)>=ni^(m)), 

m 

where /i^ff(0, 0) \ipNi) = ejvi \^pNi) and the integer N is 
the total number of particles N = n x + n y in each well. 
The index i runs from to N for each N and for the 
three lowest values of N the eigenstates are 
h/>oo) = |0),eoo = G- 

IV'lo) = jlxjOj,) ,£10 = -£• 
l^ll) = I Ox, ly) = —p» 

1^20) = |lx, lj/) ,£20 = -2£+ U xy . 

\ip2i) = ^(|2fl,,0») + l°x,2i,»,e2i = l-U xy -2fl. 

^22) = 771 (I 2 * A) - \0 x ,2 y )),e 22 = l + U xy -2ft. 

In Fig.0the mean field phase diagram has been drawn 
for the physically relevant value U xy = 1/3 which is char- 
acteristic for the proposed setup. The lobes marked MI 
correspond to incompressible Mott insulating phases with 
integer filling factors. The remaining part of the diagram, 
marked SF, corresponds to a superfluid phase with the 



columnar nematic ordering (see Fig. |SJ) discussed above. 
Considering the t = eigenstates above two things be- 
come clear. Trivially, if U xy — + the lowest lobe, and 
all other odd filling lobes, vanish and the model reduces 
to two noninteracting single flavor models as expected. 
Secondly, at U xy = 0.5 there is a level crossing between 
l^o) and |"02i) - It follows that the size of the lowest odd 
filling lobes increases with increasing values of U xy up 
until U xy = 0.5 after which it starts to decrease again. 

By considering fluctuation effects higher order in the 
tunneling amplitude, we can demonstrate that the per- 
mutational symmetry between the X- and Y- flavors can 
be broken in the Mott insulator phase. In the absence 
of tunneling, the single-particle states \ipw) and \ipn) 
are degenerate. Taking tunneling into account breaks 
this degeneracy and to second order in t (us ing for in- 
stance the Schrieffer- Wolff transformation [36j) an effec- 
tive (pseudo) spin-i Hamiltonian for the interaction be- 
tween neighboring sites can be found 

Heff = -Jeff J2 
(m.m') 

The up- and down- states of the pseudo spin operators 

(z) 

<Xm correspond to the site m being occupied by one X- 
atom or one y-atom respectively. The effective magnetic 
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FIG. 8: (Color online) Phase ordering in 3d superfluid 
phase. The directions of the arrows correspond to the phase 
angles 4> x (m y ,m z ), (f> y (m x ,m z ) and <f) z (m x ,m y ) of the or- 
der parameter fields ip x (m y ,m z ) = \4>x(rn y ,m z )\e l '^ !C ^ Tny ' m '\ 
ip y (m x ,m z ) = \ip y {m x ,m z )\e i4 ' y< ' m:c ' m ' ) and ip z (m x ,m y ) = 
\Tp y (m x ,m y )\e i ' py<m *' my) . As in the 2D (two flavors) case the 
underlying symmetry of the Hamiltonian allows for flipping 
say all the <j> x along any x-column by n to obtain another 
ground state configuration. In addition to the ground state 
degeneracy obtained from such operations, an accidental de- 
generacy associated with parallel planes of different chirality 
is present. In this figure the middle x — z plane has a different 
chirality than the other two x — z planes. 

interaction is 

_ t\Ul y + 2U xy -l) 
£ff I •: 1 E%,) 

There is a critical value of the inter flavor interaction 
U xy = V% - 1 ~ 0.414 for which J ef f vanishes. For 
Uxy > U xy , the system is ferromagnetic and sponta- 
neously favors one flavor over the other. For U xy < U£ y 
the system is anti-ferromagnetic and favors an ordering 
with X- and y-atoms on alternating sites. Thus, we 
conclude that at integer filling factor the permutational 
symmetry between X and Y flavors, (or equivalently, the 
cubic symmetry of the underlying lattice) is always bro- 
ken in the mean field ground state. Further, in the anti- 
ferromagnetic state, sublattice (i.e. translation) symme- 
try is broken as well. 

B. Phase diagram 3D, three flavors 

Using the same type of mean field theory as for the 
2D (two flavors) case, the 3D (three flavors) case can be 
treated as well. The resulting phase diagram for U xy = 
1/3 is shown in Fig. El Again, Mott-lobes with integer 
filling factors are seen surrounded by a superfluid phase 
where all order parameters ip x ,y,z have equal magnitude, 
i.e., \tp x \ = \i()y\ = \tp z \ 7^ 0. Due to the positivity of 
the coefficient U xy in the Josephson term in Eq. (|13fl the 




0.01 0.02 0.03. 0.04 0.05 0.06 0.07 



t 

FIG. 9: Mean field ground state phase diagram for the 3D 
(three flavors) Hamiltonian. 

relative phases of the three condensates are frustrated. 
Thus, writing i/j s = |-0 s |e l ^ s , s — x,y,z one finds tj> x — 

4>y = <t>y - 4>z = 4>z - <t>x = ±27r/3 ± 7r - 

An interesting effect here is that the on site frustrated 
phase configurations come in two different 'chiralities' 
that cannot be converted into each other by shifting 
any one of the phases by the 7r shift allowed by the 
Zi gauge symmetry. To see this one may consider the 
current flowing between the condensates of different fla- 
vors on a given site. The current flowing between the X 
and Y condensates on a particular site is determined by 
sin(2[0 a; — 4> y ]). In a right handed configuration with, say, 
4> x = 0, 4> y = 2n/3, <j) z — Att/3 there is an on-site current 
flowing from 

X -> Y Z X. 

The situation is different in a left handed configuration 
with cj) x = 0, 4>y = 47r/3, 4> z — 27r/3, where the current is 
now flowing in the opposite direction, i.e. 

X ^Y <- Z <- X. 

Adding an arbitrary phase of 7r (i.e. invoking the Z-i 
symmetry) to any of the phases does not affect these 
currents. 

Starting from a ground state with the same chirality 
throughout the system one can choose a set of paral- 
lel planes and change the chirality of each plane indi- 
vidually. Such changing of chirality of a plane requires 
that the whole plane has the same chirality. This ad- 
ditional ground state degeneracy is not associated with 
any symmetry of the Hamiltonian but is an accidental 
one. A similar situation occurs for special parameter val- 
ues in frustrated XY-models, where parallel zero energy 
domain walls can be inserted |2^| . One should note that 
such accidental degeneracies at the mean field level may 
be lifted by fluctuation effects associated with collective 
modes such as spin waves. 
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FIG. 10: (Color online) Example of broken permutational 
symmetry. If one increases the interspecies interaction U xy 
beyond 1/3, superfluid phases with broken permutational 
symmetry can be achieved. Shown here are the order param- 
eters sorted according to magnitude, = max(^ x ,^py,ipz), 
$3 = min(ij} x ,ipy,lj) z ) for an interspecies interaction U xy = 
0.8 as a function of t. The total number of particles is shown 
as a dashed line. 



As in the 2D (two flavors) case, the smaller Mott-lobes, 
corresponding to integer filling factors not divisible by 
the dimensionality of the system, are degenerate in the 
t = limit. This degeneracy is lifted due to tunneling, 
leading to (pseudo) magnetic ordering like that demon- 
strated for the 2D (two flavors) case. To fully lift the 
degeneracy one has to employ fourth order perturbation 
theory. The resulting Hamiltonian will include terms act- 
ing simultaneously on three and four sites. However, such 
fourth order corrections are very small and may be dif- 
ficult to observe in the proposed experimental situation. 
They can however lead to novel physics and can be in- 
tentionally generated [3H l38l] . 

Before leaving this section, we comment on the possi- 
bility of breaking the permutational symmetry among the 
flavors in the superfluid phase. As is well known, large in- 
terspecies interaction strength in the two flavor bosonic 
Hubbard model leads to phase separation. A phenom- 
ena occurring also here if U xy > 0.5. However, due to 
the positive constant in front of the 'Josephson' (flavor 
changing) term, another phenomenon can take place in 
the 3D (three flavors) model. 

As an example consider Fig. Here fi = 0.27 and 
U X y — 0.8. As can be seen for small t the system is in a 
Mott insulating state with filling factor 2. As t increases 
the system becomes superfluid. This occurs in two steps. 
First, mean field theory predicts a second order transi- 
tion to a state with only one flavor superfluid and then a 
first order transition to a state with two nonzero super- 
fluid order parameters of equal magnitude. Increasing 
the hopping strength further does not seem to make the 



third flavor superfluid. We attribute this to the large en- 
ergy cost associated with having the phases of the three 
order parameters in a frustrated configuration. 



V. INTERFERENCE PATTERNS AND 
DENSITY CORRELATIONS 

The traditional way of detecting superfluidity is by re- 
leasing the trap and looking at the density distribution of 
the expanding cloud. Provided that the cloud expands 
many times its initial diameter, the final position of a 
particle is determined by its momentum rather than its 
initial position. Hence this expanded real-space density 
distribution provides a direct picture of the momentum- 
space distribution of the trapped system. More precisely, 
the density distribution a time t after trap release is re- 
lated to the momentum density of the trapped state |$) 
as 



(n(r,t)> 



ht 



where Q(r) = mr/(ht). It is useful to think of this spa- 
tial distribution as resulting from interference of mat- 
ter waves radiated by the different lattice sites when the 
trap is released. The one-dimensional character of the 
Z 2 gauge symmetry means that thermal fluctuations can 
destroy the long range order phase order by allowing the 
phase on an arbitrary site to flip by ±7r. If the system 
disorders in this way, any interference pattern in the ra- 
diated matter waves will be destroyed as well. In this 
case, further information about the correlations in the 
system can be obtained by looking at the density fluc- 
tuations (noise) in the released cloud [3^, liol l4l| in a 
Hanbury-Brown Twiss like statistical measurement. 

We begin this section by looking at the zero tempera- 
ture momentum distribution and then consider the den- 
sity fluctuations of the expanded cloud around its mean. 



A. Interference patterns 

Although any real experiment is conducted at finite 
temperature, the zero temperature columnar phase or- 
dering may prevail for a finite system at low enough tem- 
peratures. The zero temperature momentum distribution 
is thus of interest and we will estimate it by using a single 
macroscopically occupied wave function corresponding to 
the superfluid states in the two- and three-flavor cases. 
The details of the calculations can be found in the Ap- 
pendix and we here only state the main results. 

We begin by considering a single 2D plane with N x N 
sites at zero temperature in the two-flavor system and 
model the superfluid state with a macroscopically occu- 
pied wavefunction 



(4f) 



A I 



10) 



12 



Here |0) is the vacuum state of the lattice, i.e. no atoms 
present, while a SF is the creation operator 



1 SF 



N N 

EE 



y/2N ^ ^ 

v ?n — 1 71—1 



The subscripts m and n denote the coordinates, rows and 
columns, in the lattice while a and (3 are phase factors 
(\a\ — |/3 1 = 1) determining the phase of the wavefunc- 
tion on a given site. At zero temperature the phases of 
X- particles are ordered along rows while the phases of 
Y-particles are ordered along columns, i.e., a mn = a m 
and (3 mn = (3 n (cf. Fig. and Fig. IT7|) . For a macro- 
scopic occupation M the observed density distribution 
in a single shot in the x — y plane after expansion is pro- 
portional to the momentum distribution (<i>| V'qV'Q \&) ■ 
As shown in the appendix we have for a single 2D plane 
in the two-flavor system 



<$|^qI$) = I**(Q)| 2 + |*,(Q) 



where 



= ttM 



ttM 



'o(Q) fi{Q y ,at m ) HaQx -tvk) 

odd n 

2 

'o(Q) h(Qx,Pn) Y 5(aQ y -mir) 



odd 7 



The functions 



and 



y auu </>g are the Fourier transforms 
of the on-site Wannier functions and fi(Q y ,ct m ) and 
fi(Qxi Pm) are 27r/a-periodic random functions with typ- 
ical magnitude of order unity which depend on which of 
the degenerate ground states is observed (see Fig. II 81 and 
the Appendix for details). From the above equations the 
interference pattern from a single 2D plane in the two- 
flavor system can bee seen to be a grid like structure 
as shown in Fig. ^] where the interference pattern has 
been calculated numerically for a 40x40 lattice. The ap- 
pearance of lines, rather than points as in a single flavor 
2D system, stems from the one-dimensional character of 
the superfluid state with phases only being aligned along 
rows (columns) but randomly distributed between rows 
(columns). The randomness in the distribution between 
the rows (columns) show up as the random interference 
pattern along grid lines. 

In an experiment one typically does not probe a single 
plane but it is the integrated density of a large number of 
planes that is imaged. For imaging in the plane parallel 
to the 2D planes the integrated column density (intensity 
in the absorption image) is for a N X N X N lattice with 
M atoms in each 2D plane 



I{Qx,Qy) = N 



dQ z 
2vr 



Here the line over the quantum mechanical averaging de- 
notes the averaging over the different ground state config- 
urations allowed by the Zi symmetry. Since fi(Q, a m ) = 
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FIG. 11: Calculated momentum distribution for a 40x40 lat- 
tice. The momentum distribution was calculated by numeri- 
cally summing the contributions to the distribution function. 




FIG. 12: (color online) (a) The random function 
f2(Q y ,Qz,rino) for a 40 x 40 lattice for a specific realization 
of r\n Note that fi is symmetric under inversion, (b) The 
random function g| (Qm, Qz, vino, °m) for a 40 x 40 lattice for 
a specific realization of ?j* and a m . Note that <?2 is not sym- 
metric under inversion. 



1 (see Appendix) the random interferences seen in Fig. 1111 
will be averaged out and a grid of smooth lines, void of 
interference, will be seen. Another source of smoothing 
out the random interferences comes from limited detector 
precision. For a large system, the random oscillations be- 
comes increasingly rapid and only an average over nearby 
momenta can be probed. 

In the 3D (i.e. three flavors) case, the situation is 
very similar. Special care have to be taken with acci- 
dental symmetry breaking of the ground state giving rise 
to planes of different chirality. If we assume that planes 
with uniform chirality have normals in the a;— direction 
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the momentum distribution can be written 

($1 4^q I*) = l**(Q)l 2 + l*y(Q)l 2 + l**(Q)l 5 



2Re 
2Re 



**(Q)**„(Q) 
*»(Q)**x(Q) 



2Re 



¥ W (Q)*¥,(Q) 



(18) 



where 



2n 



2 M 



f2(Q v ,Qz,Vno) 6 ( aQx 



odd ? 



i*«r = 27r| 



M 



odd n 



|* z | 2 = 27T|^| 2 y 5 |(Q a! ,Q 1 „< ) <7 TO ) ^ 

odd n 



S(aQ z — rm). 



Again, since long range order is only aligned along ID 
strips, the released cloud will be a set of intersecting per- 
pendicular planes with intersections at positions corre- 
sponding to odd momenta Qx,y,z = (2n + l)w/a. The 
planes in each direction will have a random intensity 
modulation specified by the random functions fa , g\ and 
g| (see Appendix for details). Examples of these distri- 
bution functions fa and §2 are shown in Fig. El The last 
three terms in Ea. ll8l randomlv modulate the distribution 
along the intersections of the planes. 

If a single shot measurement is made, the integrated 
column density will show a pattern of grid lines simi- 
lar to that in Fig. 1111 the grid lines showing random 
interference patterns. Between the lines a periodic ran- 
dom distribution (of lesser intensity than the lines) will 
be present. This latter distribution will be either fa or 
g% depending on the orientation of the planes with uni- 
form chirality. Thus if the absorbtion image is taken in 
the same plane as the planes with uniform chirality this 
background modulation will be symmetric under space 
inversion (cf. Fig I12f a)) whereas if it is taken perpen- 
dicular there will be no such symmetry in the random 
modulation (cf. Fig Eft))). 



B. Density-density correlations 

As pointed out above, the dimensional reduction 
present in the system means that finite temperatures 
can destroy the ID Ising-like ordering of phases along 
columns and that the individual phases at any one site 
can be flipped ±7r, i.e., the Z2 gauge symmetry is re- 
stored. In this case there will be no visible interference 
pattern although atoms are delocalized, i.e. the delta 
peaks will be smeared and a random density distribution 
will be seen each shot. 

To illustrate the usefulness of correlation measure- 
ments we consider a single N x N 2D plane in the two- 
flavor system at unit filling (M = N 2 ). If the temper- 
ature is finite, not only may the Z2 symmetry in the 



superfluid state be restored but it is also possible for the 
unit filling Mott state to be disordered. There are then 
four different possible states the system can be in 

1. Superfluid with restored Z2 symmetry. 

2. Ferromagnetic Mott insulator (all atoms of the 
same flavor). 

3. Anti-ferromagnetic Mott insulator (alternating fla- 
vors on alternating sites). 

4. Disorder Mott insulator (each site having one atom 
but with random flavor). 

If one makes multiple single shot measurements and av- 
erages the density distribution obtained in each shot, one 
obtains a measure of the average momentum distribution 
(see Appendix) 



7— i n~ r M 

(*|nq|$) = — 



I0S(Q)I 2 + I^(Q)I 5 



which is the same for each of the four states 1-4. We will 
henceforth refer to averages (■) as disorder averages. To 
distinguish the four states one can instead measure the 
HBT-like density-density correlations of the expanding 

cloud E3, E3| , 



G(r,r')E(«(rHr')) r (n(r)) t (n(r')) t 



Here (n(r)) t is the density of atoms at point r a time t 
after the trap has been switched off averaged over many 
experimental realizations (see Appendix). To measure 
(n(r)n(r')). one calculates the product of the observed 
densities n(r)n(r') in each shot and averages over several 
experimental runs. Just as for the density distribution, 
the correlation function G provides a measure of the mo- 
mentum correlations 



G(r,r') 



ht 



(nqnq/) - (n Q ) x (n Q / 



prior to trap release. 

To get a qualitative understanding of how the super- 
fluid state can be detected by correlation measurements 
we first return to the T = result in the previous section 
and look at the periodic function fa(Q y ). This random 
modulation arose because phases of Y:s were uncorre- 
cted between rows. Since the relative phases of the Y:s 
in rows is ±7r this function is even in Q and along any 
given grid line the quantity (ncjnQ/) is thus strongly cor- 
related when Q y + Q' y — 2mn/a. Averaging over many 
realizations one sees that (along a grid line of constant 
Q x — 2mt 1 'a) 



(nQTiQi) oc S(Q y + Q' y — 2mir/a) + other terms. 



In the thermally disordered superfluid state the phase 
on any site is allowed to flip by n (restoring Z 2 sym- 
metry). This destroys the delta peaks in |^(a;,y)| 2 an d 
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gives instead a random modulation given by the function 
f2(Qxi Qy)- Since only 7r-flips are allowed this modula- 
tion is symmetric f2(Qx,Qy) — fi{—Qx-,—Qy) an d we 
get strong correlations in (uqUqi) 



oc 2^<5(Q + Q' G r 



other terms. 



where G„ is a reciprocal lattice vector. On a technical 
level (see Appendix) this can be seen to arise since the 
disorder averaged propagator for single particles is nec- 
essarily short-ranged due to the random ir phase changes 
at finite temperature, while the disorder averaged prop- 
agators for pairs of particles can still be long-ranged. 

In Appendix we have calculated G(Q(r), Q'(r')) for 
the four scenarios for a single plane in the two-flavor sys- 
tem and we here give the qualitative results. The super- 
fluid state is, confirming the qualitative discussion above, 
characterized by peaks at Q ± Q' = G n (see Eq. (JEU) 
while the ferromagnetic Mott state and the disordered 
Mott state has peaks only at Q — Q' = G n (see Eqs. ijBTj) 
and (|82|) \ The correlation function in these Mott states 
are distinguished by having different background inten- 
sities between peaks and different peak strengths. The 
antiferromagnetic Mott state has peaks at half reciprocal 
lattice vectors Q - Q' = G„/2 (see Eq. (3) ). 

In the three flavor system the situation is similar to the 
two flavor scenario discussed above. In the presence of 
thermal disordering of the superfluid state the interfer- 
ence patterns of the Mott state and the superfluid states 
become indistinguishable. Again, in the three flavor case 
the pair propagators will be nonzero in the disordered 
superfluid state and G(Q(r), Q'(r')) will have peaks at 
Q± Q' = G„ (see Eq. §2$) 

While there should be no problem to measure the cor- 
relation functions for a system with three flavors, the 
two-flavor system poses a problem of technical nature 
since in an experiment several uncorrelated 2D planes 
will be created. Suppose one has N uncorrelated planes. 
If one detects one atom at position r and another at r' 
in a single experiment the atoms could have come from 
either the same plane or different planes. Measuring the 
product of densities in each shot and averaging over sev- 
eral experiments one will for the case with N 2D planes 
actually measure 



N(n(r)n(r')) t + N(N - l)(n(r)) t (n(r')) t 



rather than N(n(r)n(r')) t . Thus the signal-to-noise ratio 
scales as 1 /N requiring many experimental runs for large 
systems. 



VI. LIFETIME ESTIMATE ID 

In the previous sections, effective Hamiltonians for 
atoms in the first band(s) of the optical lattice were intro- 
duced and the mean field ground state phase diagrams 



drawn. In doing so, it was assumed that the interac- 
tion terms in the original Hamiltonian J2J) responsible 
for scattering particles between bands could be ignored. 
In this section, these interactions are taken into account 
perturbatively and the lifetime of atoms in the first band 
is estimated. The obtained (inverse) lifetime should be 
compared to other energy scales in the problem, most 
importantly the smallest one, the hopping energy. If the 
lifetime turns out to be long compared to the time scale 
of hopping, the novel states described in the previous 
sections should be possible to realize in experiment. 

To simplify matters, the discussion will be restricted to 
the ID-case. The ensuing results are expected to agree 
well, both qualitatively as well as quantitatively, with the 
2D and 3D cases to lowest order in perturbation theory. 
This follows from taking parity considerations into ac- 
count when determining the allowed transitions. Thus, 
ignoring tunneling in the y — z directions and measuring 
distance in units of the lattice spacing (see section[Hj the 
ID Hamiltonian can be written H = Hq + V with 



vox sin (0 ip n (0 



and 



V 



u 

2~ 



E 

ni ,ri 2 ,n 3 ,n 4 



d^^A^UOMO^nAO- (is) 



Here i/n(£) creates an atom at £ in the n:th band of the 
ID-system and U = AttEr i^f) Oq O (v Qv )O Q q(v 0z ). 

Apart from the field operators ipn(0 it is convenient 
to define boson operators in two other bases. First, we 
have the basis of Bloch functions with band in- 

dex n and lattice-momentum k. These functions sat- 
isfy Hou n k(0 = £ n (k)u n k(£,) and are associated with the 
field operators Snfc,a„fc- Second, we have the Wannier- 
functions <f> n (£ — 7rm) defined in section [H] In this sec- 
tion we will denote the corresponding field operators by 
a n (m) , ajj(m) . Note that these definitions depart from 
the conventions in previous sections and that operators 
corresponding to Bloch functions and Wannier functions 
are distinguished in the number of subscripts. 



A. Wide-band limit 

Begin by looking at the case when the second term 
V in eq. (|19(l is small compared to Ho and consider an 
initial state where all N atoms reside in the lowest lying 
Bloch state of the first band, 

\i) = {N\)-"\al =1 ^ /a ) |0>. 

A first order decay process is then one where two atoms 
in the first band collide, promoting one to the second 



15 



band, the other to the zeroth band, i.e. the final state is 

I ,» _ ^\i=2M^ri=0,k ^ri=l,k=TT/aO-n=l,k=TT/a ... 
U) y/N(N-l) ■ 

The first order matrix element for this transition is 



1.2 



. 8 



>=> 0.6 
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FIG. 13: (Solid line) Inverse of the coefficient w occurring in 
Eq. 12UH for the first order decay rate calculated in the wide- 
band limit. In this limit, there are no available energy states 
for the first-order decay provided 19 < vox- (Dashed line) The 
dashed line shows the maximum value of the ratio a s /a for 
which the wide-band analysis is valid. 



= U5{ko + k 2 - 2mir/a)y/N(N - 1) 

X / U 0k (0 u * 2k2 (0 u lk^/a{0 2 d^ 



If the filling factor (atoms/well) of the first band is v\, 
and the density of states of the n:th band is p(e n (k)) the 
transition rate per well becomes 



2tt { Vl Uf \J u* ko (Qu* k2 (Ou lk=7l/a (Q 2 d^\ 



p(e (ko))-l + p(e 2 (k 2 ))-i 



Defining 

w(v 0z ) = 



'S2tt 3 E r 



pM 1 +P(£2) 
this can be compactly written as 



l 0k U 2k 2 U lk=Tr/ad£, 



1 {^) 2 O 00 (v 0y ) 2 00 (v 0z ) 2 w(vo x ). (20) 



E R - sn.si 

n \ a 



In Fig. 1131 w(v) obtained from numerical calculation is 
shown. For convenience the inverse of w has been plot- 
ted. As can be seen, the lifetime goes to zero for small 
and large v. This is a result of the diverging density of 
states at the band edges. Above v w 19 the first or- 
der process is no longer energetically possible and higher 
order perturbation theory has to be applied. 



The validity of the wide-band calculation relies upon 
the assumption that the inequality v\U < t is satisfied. 
This condition can be used to obtain an upper bound on 
the ratio a s /a by assuming the lattice depth to be the 
same in all directions, i.e. vox — Vo y = v$ z = vo, which 
yields 



«s 

a 



< 



& / max 



t 



1 



This quantity is shown for filling factor v\ = 1 as the 
dashed line in Fig. IT31 



B. Narrow band limit 

From the discussion above it is clear that for deep 
enough potentials the validity of the wide-band analy- 
sis breaks down unless vi(a s /a) is extremely small. An 
alternative starting point is when i>\U ^> t while the fill- 
ing factors for the zeroth and the second band are small, 
i.e., vo, v-i <C 1. Keeping terms of order v-JJ , the relevant 
unperturbed Hamiltonian to start from is in this case 
one where tunneling events in the two lowest bands are 
completely ignored whereas interactions between atoms 
is only considered for atoms interacting with particles in 
the first band. Hence, one finds, 

h n (m) + Uox 2J no(m)n 1 (m) 

n=0,l m m 

+ ^ ra y]fti(m)(i»i(m) - 1) 

m 

+ 2uj2 I de$(0&(0&(0^»(0- ( 21 ) 

a\ l {w^a n (rn) have 



n>l 



Here the number operators h n (m) 
been introduced. 

The initial state is a product of Fock-states with defi- 
nite numbers of particles in the first band of each well as 
depicted in Fig. Here, each well m initially has n m 
atoms in the first band, i.e. 



I*> 



n v 7= io)- 



The final state is one where the population has changed 
such that, for a particular well, denoted by r, one particle 
has decayed from the first band down to the zeroth while 
another atom, in order to conserve energy, ends up in a 
Bloch state of the n:th (n > 1) band, i.e. 



l/) = 



^n r (n r - 1) 
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This state is not an exact eigenstate of the unperturbed 
Hamiltonian in Eq. l|21[l but an approximate one. The 
correction to the Bloch wave functions for n > 1, which 
will later occur in the overlap integrals, is however only 
of the order U/vo x <C 1 and can thus be ignored. What is 
more important is the associated energy shift since this 
affects the position of the band edges of the n : th band. 
This, in turn, can have impact on the lifetime since it 
affects the final density of states. Hence, one can replace 
the two last terms in the unperturbed Hamiltonian 121(1 
by a term diagonal in the band index n\ 

Ho = 51 y^£»(m)n ra (m) + U 0x ^ h (m)ni(m) 

n— 0,1 771 m 

+ ^U xx y^ y n 1 (m)(h 1 (m) - 1) 

rn 

+ ^2 ( e n(fc) + Vl\ik)n n k- (22) 

n>l,k 

Here, the first order (Hartree) shift A nk in energy due to 
interactions between an atom in the n : th Bloch band 
and the atoms in the first band have been incorporated. 

A nk = U d£ \u nk (0 + u„,- fc (£)| 2 MO 




FIG. 15: (Color online) First order lifetime it) -1 for a ID 
system with filling factor v\ = 1 and (a 3 /a) = 1/100 in the 
narrow-band limit according to Eq. 112^1 . The solid lines show 
the ratio between the lifetime and the time scale for hopping 
(M^ 1 ) for n r particles in well r. From top to bottom; (blue) 
n r — 2; (green) n r — 3; (red) n r — 4; (black) n r = 5. The 
dashed line shows the ratio ti/v-JJ which should be less than 
unity for perturbation theory to be valid. The dot-dashed 
line shows the result obtain by using the wide-band formula 
in Eq. I|2()|l using the same parameters. 




FIG. 14: Typical initial state \i) for the lifetime estimate 
in the narrow band limit. All atoms are residing in the first 
band, localized in the wells. In this particular case the filling 
factors are Vq = 0, v\ = 1 and v n >\ = 0. 

For first order decay one needs only the matrix element 
(/| V \i) from which the rate follows; 



2tt 2 , 
— U n r (n r - 



1) J d£Kk(0m)M0 2 P{tn{k)). 

To present a comprehensive numerical analysis of this 
decay rate is prohibitive due to the large number of pa- 
rameters entering expression Eq. J23J). Thus, for sake of 
illustration, we will here restrict the discussion to unit 
filling factor in the second band, i.e. v\ — 1. Further, 
we use (a s /a) = 1/100 which is a reasonable value from 
an experimental point of view. The lattice depths in the 
transverse directions will be chosen slightly larger than 
in the x-direction, i.e. chose vq v — vq z = vq x + 1. 

The results of the calculation, using wave functions 
obtained from band-structure calculations, are shown in 
Fig. El which plots the ratio between the hopping rate 
and decay rate, t\/{%w). The different solid lines cor- 
respond to different number of particles initially in the 
well. The cases n r = 2,3,4,5 are shown, n r = 2 hav- 
ing the longest lifetime and n r — 5 having the shortest. 



The dashed line shows the ratio t\/(v\U) which should 
be less than unity for the expression to be valid. As a 
comparison, the resulting lifetime obtained in the wide- 
band limit Eq. i12()ll is also shown as the dash-dotted line. 

The most interesting part of the result shown in 
Fig. (|15fl is the sudden decay of the lifetime. This is, 
as was the case in the wide-band limit, a result of the 
diverging density of states pfaik)) near the band edge. 
For lattice potentials deeper than vo x w 20, there is no 
phase space (no available final energy levels for the ex- 
cited particle) , available for the first order decay. To find 
out the lifetime for larger values of vq x , second order per- 
turbation theory is needed. 

Consider again the same initial state \i) as above. 
Adhering to energy conservation arguments, there are 
three different, mutually orthogonal, final states reach- 
able through a second order process, 



I. A 



«ifc«o( r ) 2 "i( r ) 3 

y/2n r (n r — l)(n r — 2) 



l/ 2 ) 



1/3 



&L4(r) 3 aiM 4 
y^6n r (n r — l)(n r — 2)(n r — 3) 



a i'fc' a ifc a o(r) 2 Qi(o 4 

yj2n r (n r — l)(n r — 2)(n r — 3) 
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FIG. 16: (Color online) Second order lifetime w^ 1 for a ID 
system with filling factor V\ = 1 and (a s /a) = 1/100 in the 
narrow-band limit according to Eq. 123H . The solid lines show 
the ratio between the lifetime and the time scale for hopping 
(M^ 1 ) for n r particles in well r. From top to bottom; (green) 
n r — 3; (red) n r = 4. 



The corresponding decay rates ^1,2,3 are obtained from 
the text-book relation 



2tt 



y, (fi\V\m) (m\V\i) 



Numerical evaluations of the decay rates reveal that the 
dominant contribution to the total decay rate Wtot = 
W1+W2 + W3 comes from w\ . The reason for this is easily 
understood; the contribution from wi is small due to de- 
structive interference of time reversed processes while the 
smallness of w$, is due to smallness of overlap integrals, 
which in turn can be understood from parity considera- 
tions. 

The decay rate wi is shown in Fig. ^| As can be 
seen, it is possible to achieve lifetimes considerably larger 
than the inverse of the hopping energy, thus justifying the 
validity of the Hamiltonians in eqs. (|11|) and i|13|) . 



VII. CONCLUSIONS 

By extending the usual mapping to the bosonic Hub- 
bard model of ultra cold atoms in an optical lattice 
to incorporate higher Bloch bands, effective Hamiltoni- 
ans governing the dynamics of atoms in the first Bloch 
band(s) have been obtained. These Hamiltonians resem- 
ble previously studied bosonic Hubbard Hamiltonian but 
differ in two important respects: 

• Atoms in the first excited band are labelled by three 
possible flavors X, Y, Z. The dynamics is such that 
X particles can (to a good approximation) move 
only in the x direction, etc. 



• Flavor changing collisions of atoms on the same site 
leading to conversion of the form XX — > YY, 
etc. occur. 

By appropriate choices of the lattice depths in the dif- 
ferent directions the number of flavors and the effective 
dimensionality (equal to the number of flavors) of the 
system can be changed. To obtain values of the rele- 
vant parameters, such as hopping energy and interac- 
tion energies, entering these effective Hamiltonians we 
have solved the time independent Schrodinger equation 
(Mathieu equation). 

The effective Hamiltonians in two and three dimen- 
sions also show, apart from the usual global (7(1) gauge 
symmetry, a set of Zi -gauge symmetries intermediate be- 
tween local and global. The ground state in the 3D (three 
flavors) case also displays a chiral symmetry breaking and 
an additional accidental ground state degeneracy associ- 
ated with different planar chiral ordering. 

The phase diagrams for two particular cases relevant 
for experiment have been sketched using mean field the- 
ory, indicating quantum phase transitions between Mott- 
insulating and superfluid states. 

Using time dependent perturbation theory up to sec- 
ond order in the interatomic interactions the lifetime of 
the atoms in the excited bands have been estimated. 
The results show that life times considerably longer (or- 
ders of magnitude) than relevant dynamical time scales 
can obtain. This suggests that it may be possible to 
realize quasi-equilibrium in the subspaces of meta sta- 
ble states spanned by the effective Hamiltonians. Fi- 
nally, we would like to stress that the mean field the- 
ory used to draw the phase diagram is only able to de- 
scribe the most simple scenario with a transition from 
a Mott-state to a superfluid state with order parame- 
ter (X) ^ 0. It is well known H E E H El that 
other multi-flavor bosonic Hubbard models such as the 2- 
species Bose-Hubbard model shows a rich phase diagram 
with phases which cannot be described in this simple 
approximation. The present model, already rich at the 
mean-field level warrants further study. In particular, we 
have pointed out potential connections to certain classes 
of models of frustrated spins 0, |U 0] an d bose 
metals j^] that also have an infinite but subextensive 
number of Z2 gauge symmetries and as a result exhibit 
dimensional reduction and exotic phases. With the mi- 
croscopic Hamiltonian developed here, these connections 
can and should now be pursued in detail. 
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APPENDIX 

Here we provide a detailed derivation of the density 
distribution one expects to observe in the various phases 
and the different ways of measuring the density density 
correlations in the released cloud of atoms. 

If the system is in a many body quantum state |$) 
when the trap is released at time t = the density dis- 
tribution of atoms at a later time t is given by 

(n(r)> ( = ($|[/t(t)n(r)[/(i)|$), (24) 

where U(t) is the time-evolution operator of the released 
system U(t) — exp{—ihT 1 Ht). To measure the quantity 
in Eq. I|24l) one has to, in general, perform several mea- 
surements starting with the same trapped state | $) each 
time. An exception to this is when the ground state is to 
a good approximation a macroscopically occupied single 
particle state. This is typically the case for a superfluid 
system and a single measurement gives a good approx- 
imation of (n(r)) t . For a weakly interacting dilute gas 
of atoms the interactions between atoms can be ignored 
during the expansion of the cloud and the time evolution 
operator in Eq. (|24|l can be replaced by the free time- 
evolution operator Uo(t). Expanding in the momentum 
components one finds 

f n ( r )\ - [ dkl f rfka r -»(k 1 -k 2 V(r-^(k 1 +k 2 » 

x ($|4^ k2 |$). (25) 

For a system of linear size L and for times ht ^> ml? the 
stationary phase approximation gives 

(HI \ 3 TTI V 

-) ($|n Q(r) |$) Q(r) = ^. (26) 

Measuring the density of atoms after a long time of flight 
t thus corresponds to a measurement of momentum dis- 
tribution of the state |$) prior to trap release. 



In a typical experiment one takes an absorbtion image 
of the released cloud. This means that the only the Inte- 
grated column density is measured, i.e., if an image of, 
say, the x — y plane is taken, one measures 

I(x,y) = J dz(n(r)) t = (g) 2 J ^ <$| n Q(r) |*> . 

In the next subsection we derive the momentum distri- 
bution ($| 7"iq( r ) |$) for the superfluid states in the two- 
and three-flavor systems at zero temperature where Z 2 
symmetry is broken. 



1. 2D, two flavors, superfluid state, T=0 

For the two-flavor case the system is comprised of 2D 
planes with uncorrelated ground states. A superfluid 
state of a single 2D plane can be described by a wave 
function with M particles in a single state 

|$) = (M!)- 1 /2( a ^ |o>. (27) 

^ N N 

a SF = /q at- { amn ^-rnn + PrnnXmn) • 

m — 1 n— 1 

The subscripts m and n denote the coordinates, rows and 
columns, in the lattice while a and /3 are phase factors 
\ a \ = \P I = 1 determining the phase of the wavefunction 
on a given site. 

To evaluate the momentum distribution we expand the 
field operators and ipQ m terms of the localized cre- 
ation and destruction operators X^ nn ,Y mn etc. where 
the subscripts m and n respectively denote the row and 
column for the site on which the operator is acting. For 
a general state |$) (not necessarily the state in Eq. I|27|)) 
we find 



iQ-(ri— r 2 ) 



x ($| [xl^Xnn^iT +Yi in y mini {v 1 )*] [X m2n2 <p* m2n2 (r 2 ) +r m2 „ 2 e 2 „ 2 (r 2 )] |$) . (28) 
The localized Wannier orbitals <pnm ( r ) an d 't'nm ( r ) can be rewritten 



CM = (-l)'VS(r - nax may), ft„(r) = (-l) m $(r 



nax — ma 



y) 



with the prefactors (— l) n ( m ) coming from the gauge choice in the inital way of writing the Hamiltonian in equation|21 
Carrying out the Fourier integrals we find 



$i^ q i*> = E E e iQ(Ri - R2 > (* fe ini (-i) ni ^(Qr+^ ini (-i) mi ^(Q) 



mini 7ri2n2 



[A m2 „ 2 (-l)"^g(Q)+r m2 „ 2 (-l)" l2 ^(Q) 



4> 



(29) 
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Here the position vectors Ri and R2 are shorthand for 
the lattice vectors 

Ri = n\ax + miay, R2 = n^ax + m-2.ay 

and <^q^(Q) denote the Fourier transform of the onsite 
wavefunctions. In the Harmonic oscillator approximation 
these are given by 



_ _3/4 -5/2, 



6g(Q) = ^/ 4 7 



t2 i 1 ) 2 i 1 ) 2 



2ir\/mVo 



nx 



(30) 



To evaluate ($| V'qV'q |^) we need to calculate the expec- 
tation values of the kind ($| X^Y |$). For the superfluid 
state |$) in Eq. (|27|) it is easily verified that in terms of 
the in-plane density p = M/N 2 one gets 



Y 



^) — 2 a ™ini am 2n2 



efc. Hence the terms in Eq. <|29[1 factors and one can 
write it conveniently as 



(#|^qI*> = l^(Q)l 2 + l*y(Q)l 5 

+ 2Re [f x(Q)**y(Q) 



(31) 



where we have defined 

<MQ) = ^(Q)7|^ e -^- R -(-l)"a mn (32) 

mn 

*,(Q) = ^(Q)J|^ e - lQ - R -(-l) m /3 mn (33) 

mn 

For a system at absolute zero the phase factors a and 
/3 are aligned along rows and columns respectively but 
are randomly distributed between the lines and columns. 
To describe this situation we introduce two sets of fields, 
77^ and 1]% which can take on values ±1. The relation 
between these values of the fields and the phases along 
rows and columns is shown in Fig. 1171 Thus we can write 

Consider now the summations needed to evaluate ^ x 



*x(Q) = ^(Q)^E e " iQ ' Rm "(- 1 )"^ ( 35 ) 



The summation over columns (n-summation) converges 
in the large TV limit to a sequence of delta-functions 



(36) 



»1 = 1 



€ 1=1 



■''n+2 = 



>4 = i 
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FIG. 17: Sample configuration of phases and the fields rf^ and 
77^ for a plane in the two-flavor system at zero temperature. 



and a similar equation can be obtained for ty y . Hence 



= 2ttN!- 



5 (aQ x - (2n + 1» 



v „-«(m-m')0,o n i 2: 

^ / j c 'Im'lm' 

mm' 



Introducing A = m — m! the last summations can be 
rewritten 



Nfl(Qy,Vm) 



\^ -i{m-m')Q v a x x 
° 'Im'lm' 



mm' 



E e " AQya E«-A 



N + E e- lAQ » Q ^ 

A/0 m 



^m-A (37) 



where we have defined the random momentum distribu- 
tion function /1 . With the aid of Eq. lf3*?|> we now deduce 
some properties of f\. We begin with the magnitude of 
the function for any value of Q y . For each nonzero value 
of A the last summation is over an uncorrelated sequence 
of integers ±1 and can be viewed as a ID random walk 
for which we have that 



0(VN) 



The summation over A contains N — 1 terms which for 
each value of Q y are random of magnitude \/~N. This is 
again a random walk with N — 1 steps and we conclude 
that the whole expression in Eq. I|37|) is of order N. This 
can also be seen by noting that 



f [ n/a dQ y f 1 (Q y ,r,l) = l. 



Thus for each configuration rf^ we have a randomly oscil- 
lating function f(Q y , r] m ) of unit magnitude. An example 
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FIG. 18: Example of the random function fi(Q y ,rim) defined 
in Eq. . 



The subscript j denotes collectively the x- y- and z- 
coordinates in the 3D lattice. As in the two flavor case, 
the observed momentum distribution can be written as 



with 



l**(Q)l a + l*„(Q)l 

2Re [* X (Q)** 3/ (Q) 
2Re [* y (Q)** z (Q) 
2Re [* Z (Q)**,(Q) 



l<MQ)l 5 



(38) 



**(Q)MS(Q)J?£' 



,-iQ R„ 



of fx obtained for a specific realization of 77^ with N = 40 
is shown in Fig. ED From Eq. (|3T|) it is also clear, since 
77^ and r? m -A are uncorrelated for nonzero A that the 
average over allowed ground state configurations is 



h(Qy,r^) = ^ £ / 1 (Q„0 = 1. 

»?i '•••>')n= ±1 

An important property of f± is that it is even in Q y 

This is a result of the nematic ordering between rows. 
I^yl 2 can be calculated the same way as \^ x \ 2 and we 
get 



odd n 



odd m 

The interference term, the last part of Eq. (|31|l . for the 
momentum distribution vanishes. To see this one can 
make use of equations <|5I|) and 



Vx 


2 

= ttM 






2 


0o (Q) 




= TTM 



2Rc 



**(Q)**«(Q) 



= 47^00 0o Re 



X z 



mn 

<5 (aQ y - (2m + l)7r) e" 



'In 



2. 3D, three flavors, superfluid state, T=0 

For the three flavor case at T=0 we consider again a 
state of the kind in Eq. (|27|l but with 



L SF 



N 3 



*„(Q) = $(Q), 



mno 



* 2 (Q)^(Q)^£ e -Q' R " 



'(-i) n A, 



-1)°7* 



Here the subscripts mno refer to the y— and z— 
coordinates in the lattice respectively. To see how to 
handle the phase factors in the three flavor case, we begin 
with a state without accidentally broken chiral symmetry 

O^mno — Vnoi Pmno — e ^ Vrnoi Imno — e Vmn 

where the random fields rjij can again take on values ±1. 
Since the accidental chiral symmetry breaking occurs in 
parallel planes we can without loss of generality single out 
the x-dircction as the direction in which planes have uni- 
form chirality (To compare with Fig [S] make the rotation 
of axes y — > z, z — > x, x — > y in Fig. |SJ. We thus intro- 
duce an additional field a m taking values ±1 for planes 
with different ir-coordinate m. The corresponding phase 
factors for such a state will be 

Oimno = Vnoi Pmno = £ 3 Vmo> Imno = B m 3 T] mn (3Q) 



We can now evaluate 
two flavor case 

2 - 2 M 

T 



in the same way as for the 



2tt 



o 



f2(Q y ,Qz,Vno) £ 5 (aQ x - mir) 

odd m 

(40) 

In Eq. lj4Tj|l f2(Q v ,Qz,Vno) nas been introduced 



-n 2 )Q y p ia(o 1 -o 2 )Qz x x 

'm 1 oi'm202' 



n 2 o 2 



The random distribution function / 2 is the two-variable 
analog of the function fi above. An example of f% for a 
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40x40 lattice is shown in Fig. l^f a'). Just as fi, fa obeys 
a sum rule 



(— 

V27T 



7r / a pTr/a 



dQ y dQ z f 2 (Q y ,Qz,Vno) = 1 



-7r/a J — n/a 

is symmetric under inversion 

f2(Qy,Qz,VZo) = h(-Qy,-Qz,Vno) 

and has an average equal to unity when averaged over 
ground states 



h(Q y ,Q z ,vlo) = l- 

2 



The expressions for 





2 






and 





are similar but the 



accidental ground state degeneracy modifies the random 
distribution functions. Explicitly we have 

|*„| 2 = 2n\4>l\ 2 ^-g y 2 (Q z ,Q x ,7 1 y m ,a m )J2S(aQ y -nn) 

odd n 



\^ z \ 2 = 2TT\4> z \ 2 ^-g^(Q x ,Q y ,^ n ,a m ) ^ S(aQ z - mr) 

odd n 



with 



n V = „ja[(mi-m 2 )Qx + (oi-o 2 )Q i ] y y .-ifK -"m S ) 

i/2 — / y c '/miOi'/m 2 o 2 c 



.'72 



E io|(rai-m 2 )Q I + (ni-ri2)Q,] z z - 

° '/mini '/m 2 n 2 ° 



: (tT m i-0- m2 ) 



mini 

m 2 n 2 



En example of the distribution function g| is shown in 
Fig. I12f b) . Note that due to the fields a characterizing 
the different chirality of planes this distribution function 
is not symmetric under inversion. Finally we look at the 
interference terms in Eq. i|38|) . 



** x * y = 47r 2 0g^| ]T «(oQ«-m7r)*(oQ B -n7r)X; E e ^" C ^ ) (- 1 ) ni <o 1 (-l)^< l0a e i * <rma 

rn.n odd ri\0\ m^o-z 

(41) 



In the above equation the summations over n\ and mi 
constitute random walks. For the ri\ summation this is a 
random walk on a line with N unit steps ±1 giving rise 
to, for each value of 0\ a random term of order y/~N. The 
sum over 777,2 can also be viewed as a random walk for 
each value of 02 but in the complex plane. Each step be- 
ing of unit length in any of the four directions ±27r/3 and 
±47r/3. Summing over ri\ and 771,2 thus yields, for each 
(01, 02) a random term of magnitude N with a completely 
random phase. Thus the interference terms in Eq. I|38() 
(the other two terms can be treated similarly) give rise 
to a three dimensional grid of lines in the released cloud 
where the density along any given line is randomly dis- 
tributed. If the density is averaged over several shots, 
with different ground states we have no contribution from 

the interference terms since ^ x ^ y = 0. 



3. Density averages and correlations , T > 

If T is large enough for thermal fluctuations to restore 
Z2 symmetry but still small enough to preserve the dis- 
tinction between the Mott state and the superfluid state, 
measuring the density distribution alone does not suffice 
since the delta-peaks will be smeared. Instead correla- 
tions can be measured. To this end, assume we have a 
single physical system. At finite T this system undergoes 



transitions in a manifold of M states. Denote this man- 
ifold by the states In a single shot a single 
one of these states will be probed. In an infinite series 
of experiments each of these states will be probed an in- 
finite number of times and one can thereby measure the 
quantity 



($1 6m J2 mom. m 

i—l 

Here we have ignored the Boltzmann factors since the 
manifold we are looking at is nearly degenerate. In reality 
only a finite sequence of M. experiments can be carried 

out and the fluctuations in ($| O is of concern. There 
are two sources fluctuations; First, for each state 
there is quantum shot noise. Second, since not all of the 
N states will be probed there will be deviations due to 
not sampling the entire distribution. 

If the manifold of states probed are superfluid states 
then Af = 0{[2d] Nd ) with d being the dimensionality of 
the system. Since superfluid states are to a good approx- 
imation macroscopically occupied single particle states, 
fluctuations due to shot noise are reduced. The remaining 
fluctuations are classical and expected to scale as Mr 1 ! 2 

and ($| O |$) should in principle be possible to measure. 
On the other hand, if the state measured is a Mott state 
the manifold consists typically of only a few 
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states in which case multiple measurements reduces the 
quantum shot noise since each quantum state will be 
probed many times. We thus conclude, that by making 
repeated measurements and averaging the results, one 

can measure (3>| O |<&). 

A quantity of interest to measure in this way is the 
correlation function 



G(r,r') = (n(r)n(r')) i - (n(r)) t x (n(r')) t . 

Again, if %t ^> ml 2 this is to a good approximation the 
same as 



G(r,r') = ( 



ray 
hi) 



(n-QUQ,) - (nq) x (tiq/) 



The disorder-averages of the momentum density distri- 
butions are easy to calculate. For instance, for the two 
flavor superfluid state in Eq. (|27|l we have 



l**( Q)l 2 + hMQ)l 2 
2Re [* a! (Q)**„(Q) 



(43) 



The averages in Eq. (|43[) can be calculated using the rep- 
resentation in Eqs. I(32|l and 1)33(1 



I*, 



The average means averaging over all a mn = ±1 (and all 
(3 mn = ±i). It follows 



M, - 



l^(Q)l 2 = yI0S(Q)I 2 l**(Q)l 2 = T l0o y (Q)l 2 



and ^ r x (Q)*^j / (Q) = 0. Hence, for the two flavor case we 
find 



(*I«Q|*> = y(l^(Q)| 2 + l0o y (Q)l 2 ), (44) 
whereas for three flavors we have 

W^W) = y (l^S(Q)l 2 + l^o(Q)l 2 + I0S(Q)I 2 



We now turn to the evaluation of the two point corre- 
lator which we begin by normal ordering 



i Q n Q ,) = (2n) 3 (n Q )S(Q - Q') + U^Mq 



The normal ordered expectation value ^V'qV'q'^Q^'Q' 
can be written in a form analogous to Eq. (|29|l 



e iQ-(Ri-R3) e iQ'-(R fc -Ri 



ijkl 
X /$ 



X 



[x](-l)»^g(Q)*+F i t(-l)^^(Q)*] [A^(-1)"^5(Q')* +n t (-l)" ^fc ^(Q') , 

x i (-i)"^5(Q) + y J -(-i)^^(Q)l [^(-^"'^(QO + ^C-ir^CQOl I* 



(45) 



r 



Here the subscripts ijkl are collective row and column 
coordinates for the site index in the 2D lattice. 

For the two flavor superfluid state in Eq. I|27|l (a single 
plane with N x N sites having a total of M particles 
it is easy to verify that the expectation values of on-site 
operators are given by expressions of the type 

To calculate the average over disorder we have to average 
over ctj — ±1 and (3j — ±i. The nonzero averages are 
easily seen to be 

(46) 



81 

4 



= -j [SikSji + SijSki + Su5 kj ] (47) 



A -a*alf3 J f3 l = - P -k^i (4si 



nn^i) = ^ a * kPjai = 4^ 4; (49) 



) (YfxlXtf) = P-Plala^^^-daSk, (50) 

XjY^Xt) = ^atPtfrat = £t a 5 kj (51) 



(*M x i Y i) = x a ^ a ^ = T 4 '^' (52) 



yMxjXi) = Y^^ ai = -T 6ik5jl (53) 



[SikSji + SijSki + SuSkj] ■ (54) 



Note that Eqs. l|i5 )l .l|i7 )l .l(55 )l and JS!} have contributions 
that correspond to pairs of particles propagating. The 
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disorder average of single particle propagation being zero 
due to the random orientation of phases. Using Eqs. 145|) - 
H54JI we can evaluate the terms in Eq. (|45(l which are 
nonzero. We state each term contributing to the correla- 



tor in Eq. (|45|l separately using subscripts to denote the 
specific ordered combination of operators from which the 
term derives. 



$i^q^qv<q'I<i>; 



xxxx 



(4> 



V'qV'q^qV'q' 



' YXYX 



4>) 



'yrxx 



yyyy 



6g(Q)| 2 |<^(Q')| 2 ^ £ [1 + eiCQ+Q'MRi-Rj) + e <(Q-Q').(R«-it, 

•^(Q)>o I (QT^(Q)^(Q')jE el(Q+Q ' HR " Rj) 



4 
M 2 

*S(Q)>S(Q')>S(Q)^o(Q')— 

2 



<$>) YXXY = 0o s (Q)>5(Q')>§(Q)^(Q')xE ei(Q_Q ' HR ^ R 



^(Q)*^(Q')>[!(Q)0S(Q')^E ei(Q " Q ' HR ^ R 



M 2 

= ^(Q)>g(Q')>o(Q)^(Q')-r 



6g(Q)>g(Q')*0g(Q)0S(Q')x E e* (Q+Q ' HR *- R 



W)I 2 I^(Q')I 2 ^E[ 1 + 



e i(Q+Q')-(Ri-Rj) _i_ e i(Q-Q')-(R<-Rj 



(55) 

(56) 

(57) 
(58) 

(59) 

(60) 
(61) 

(62) 



Collecting the results of Eqs. ll4^t- ll^|) we find 

,M 



G 2 s £(r,r') ex (2^) 3 T ^(Q-QO(|^(Q)| 2 + |^(Q)| 2 )+^|^(Q)0g(QO + 0g(QK(Q')| ^e^^'H^ 



R, 



4 



KQ)^S(Q') " $(Q)$(Q')| ^eKQ+QO-CR.-R. 



i) 



(63) 



where the factor of proportionality is (S) ■ The Fourier 
sums give, in the limit of an infinite lattice, sequences of 
delta functions 

£ e i(Q±Q')-(R i -R 3 ) ^ f 2 £ jqq ± q/] _ G .) 

ij i 

where are reciprocal lattice vectors. The most inter- 
esting part of Eq. (|63|l is the second line which comes 
from the pair like propagation. This can be used as a 
signature to detect the superfluid phase even if thermal 
disorder has restored the symmetry. 

For comparison we also look at the 2D (two flavors) 
Mott state. For simplicity we consider unit filling. There 
are three scenarios for the unit filling Mott state that 
need to be considered; (a) FM Mott state, i.e. all atoms 



of the same flavor (b) AFM Mott state, X-flavor and Y- 
flavor on alternating sites and (c) thermally disordered 
Mott state with random occupation of X- and Y- flavor 
on each site. 

In the ferromagnetic Mott state at unit filling M — N 2 
we have two degenerate ground states |$i) = J\ i Xj |0) 

and |$ 2 > = Ui Y i 1°) and the average in Eq. j|4"2*|) is trivial 
to evaluate 

($|n Q |$) = -(($i|n Q |$i) + ($2|n Q |$2)) 

= y(l^(Q)| 2 + l^(Q)| 2 )- (64) 

The momentum correlator can be calculated using 
Eq. I|45|) . For the state this equation reduces to 
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($i|^,^Q^Q'l*i> = |^(Q)| 2 |^(Q , )| 2 Ee^^- R i)e i Q'-^- R ')($ 1 |XtxtX J .X J |$ 1 ). (65) 

ijkl 

There are two pairings of operators that contribute to the average 

(^ixjxlXjXt |*o - (i - * ifc )(Mw + <^M- 

The term (1 — 8ik) results from having only one particle at each site but we will ignore this term (and terms similar 
to it in what follows) since its relative contribution is of order 1/N 2 . The disorder average contains only two states 
yielding 



<*l ^UqV>q' I*) = ~ (|05(Q)I 2 I^(Q')I 2 + |0 Q y (Q)l 2 l^(Q')l 2 ) £[i + eW-Q'MR,-**) 

ik 



(66) 



The correlation function for the Ferromagnetic Mott state is thus 



M M 2 
G% D M (r, r') « (2^) 3 <5(Q - Q') y (|0§(Q)| 2 + l^(Q)| 2 ) + — (I0§(Q)| 2 - l^(Q)| 2 ) (I0§(Q')| 2 - l^(Q')| 2 ) 



I fe(Q)l 2 l^(Q')l 2 + l0o y (Q)H^(Q')l 2 ) E e ' (Q_Q ' HR " R 



(67) 



In the antiferromagetic Mott state the disorder aver- 
age is again over two states. Dividing the 2D lattice 
into two sublattices A and B these states are |$i) = 

n ig An i6B *Wio> and 1*2) = n ie xn igB ^io>. 

For the momentum density we have 

Mr- 

(*iKI<fi) = (* 2 |n Q |$ 2 ) = T [|0g(Q)| 2 + |^(Q)| 2 



hence we find again that 
M 



($|nq|$) 



l<^(Q)l 2 



TO)I 2 



The normal ordered two point correlator can again be 
written in the form of Eq. I|45|l and has 6 nonzero contri- 
butions. The disorder average over the two states will in 
this case make no difference since the two different states 
always give the same contribution and it is enough to 
consider one of them. 



($i|^^Q^Q'I$i)x,h = l^(Q)| 2 |^(Q')| 2 (^-+ E < U g " " 

\ ieAjeA 

<*i| ^ Q *Mq>\*i )yyyy = W(Q)| 2 |^(Q0! 2 ( ^ + E e^-Q'H^. 



-it,) 



ieBjeB 



XYXY 
XYYX 

(^iIV'qV'q^qV'q' |$l) YXYX 

YXXY 



M 2 



l$(Q)l 2 l$(Q' Ma 



0g(Q)*^(Q')*^(Q)^(Q') E eW-<*>^- 

ieAjeB 

M 2 ~ 

— |^(Q)| 2 |0S(Q')| 2 



R, 



$(Qr<£§(Q')>o(Q)$(Q') E e i(Q - Q ' HR '- R< 

Hence, we find for the 2D antiferromagnetic Mott state at unit filling the correlation function 

Mr- 

G 2 A D FM (r,r') « (2^) 3 ,5(Q-Q') T [|0S(Q)| 2 + |^(Q)| 2 



(68) 

(69) 

(70) 
(71) 

(72) 
(73) 
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+ 2Rc [^(Q)Vg(Q')Vg(Q)^(Q')J E cos t(Q - Q') ' ( R * - R ^ 

+ (|0S(Q)| 2 |^(Q')l 2 + l^(Q)H^(Q')l 2 ) E eW-Q')-^-^). 



iEAjEA 



The Fourier sums converges in the limit of large iV to 



E cos [(Q - Q') • (Rj - Rj)] = ^ EK- 1 )" + ~ V )5(Qs " T } 



ieAjeB 



(74) 



E « <(Q - Q ' ,,(B '- R,) = ^ Et 1 + (-IJ^W- - ™ V } 

ieAjeA mn 

and the correlation function for the antifcrromagnetic Mott state will thus have peaks at locations corresponding 
to half reciprocal lattice vectors. ^ 



Finally we look at the disordered Mott state where each 
site holds one atom but whether it is an X or a Y is 
random. The manifold of states to average over thus 
contains Af = 2 N states. Such a state can be written as 

\*) =Ul[x}(l + Vi) + YHl-Vi)} 10) 



where r\i is a random field taking on values ±1 on each 
site i. The disorder averaged momentum distribution is 



again the same as before 



<d>i4^Qi$> = E^^i(Q)i 2 + 



' f " :J ' ^(Q)l 2 



M 

Y 



(^l(Q)l 2 + l^(Q)l 2 ) (75) 



and there are six contributions to momentum correlator 



m^ Q 'Mc,m xxxx = ^s(Q)i 2 i0S(Q')i 2 E 1 + et(Q-Q ' HR ^ Rfc) 



ik 



WW q >Mq>\*) yyyy = Jl^(Q)l 2 l^(Q')l 2 E 1 + e4(Q_Q ' HR ^ Rfc) 



($| Vq^q'V'qV'q' I*) 



XYXY 



M 2 - 

— |^(Q)| 2 |^(Q')| 2 
1 - 



XYYX 



§(Q)*^(Q')*^(Q)^o(Q / ) E e^Q-Q')'^-^) 



($| Vq^q'V'qV'q' I*) 



YXYX 



4 

M 2 - 

— I^(Q)I 2 I0S(Q')I 2 

1 7 



ik 



YXXY 



T^(Q)>§(Q')^S(Q)^(Q')E el(Q_Q ' HR ^ R 



■*) 



ik 



resulting in a correlation function for the disordered Mott state 



G 2 D D (r,r>) cx (2^) 3 «5(Q - Q')y [|<^(Q)| 2 + l$(Q)| 2 



1,7 

4 



(76) 
(77) 

(78) 
(79) 

(80) 
(81) 



l^(Q)^5(Q') + ^(Q)*^(Q')l 2 E eJ(Q_Q ' HR ^ Rfc) - 

(82) 



ik 



4. Correlations 3D, three flavors, T > 



uate the correlation function (<&| V'qV'q'^qV'Q' 1^) onc 



We now look at the momentum correlations in the ther- 
mally disordered three flavor superfluid phase. To eval- 



can write down the extension of Eq. Q45[l. There will 
be overall 81 terms in the expansion to evaluate. When 
taking the disorder average only 21 terms are nonzero. 
Note that when taking the disorder average in the three 
flavor model one has to average not only over all possible 
7r flips of the phases but also over the symmetry breaking 
field a m (see Eq. to account for the chiral symme- 

try breaking as well as over the 3 directions directions in 
which chiral symmetry is broken. The nonzero averages 
one obtains in this way are shown below 

( x l x l X 3 X ) = ( Y i Y k Y 3 Y l) = ( z l z l z i z i) 

2 2 

= -Tri^ikSji + Sij5 k i + SuSkj] (83) 
( x l x l Y j Y i) = (Xlxlz^t) = —^SfkSn (84) 

( y M x j x i) = ( y M z j z i) = -\Y 5ikSjl (85) 
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= ( z 14 x j x l) 


) = 


(86) 


x M x j y i) 


= ( x \ z k x i z l ! 


> = ^ 5 ij5ki 


(87) 


Y ? x l Y j x i) 


= ( Y }4 Y i z l) 


= -^Sijhi 


(88) 




= ( z M z j y i) 


= -^kjhi 


(89) 




= ( x I4 z j x i) 


) = ^-gdilfikj 


(90) 




= ( y I4 z j y i) 




(91) 


z l x l x j z i) 


= ( z M Y i z i) 


= -y5u6 kj 


(92) 



and the desired correlator is obtained 



G 3 s ^(r,r') ex (2^) 3 A(Q-Q')^ [|^(Q)| 2 + |^(Q)| 2 + |^(Q)| 2 ] + ^ (l^(Q)^(Q') - $(Q)$(Q') 
+ |^(Q)^(Q') - 0S(Q)0g(Q / )l 2 + l^(Q)^(Q0 - 0g(Q)^(Q')| 2 ) E e j(Q+Q ' HR '- Rfc) 



9 



|^(Q)^(Q')* +^(Q)^(Q0* +^(Q)^(Q')1 2 E el(Q_Q ' HR ^ Rfc) - 

I 



(93) 



ik 
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